#ifndef UHM_MAT_EXTENSION_LU_PIV_HXX_
#define UHM_MAT_EXTENSION_LU_PIV_HXX_

namespace uhm {
  
  // Basic LU piv...
  template<class A_, class B_>
  class LU_piv_ : public LU_nopiv_<A_,B_> {
  private:
  protected:
  public:

    LU_piv_() : LU_nopiv_<A_,B_>() { }
    LU_piv_(Int_ type, Int_ uplo, Int_ sym,
            Int_ mt, Int_ mb,
            Int_ rhs, Int_ blk) 
      : LU_nopiv_<A_,B_>(type, uplo, sym, mt, mb, rhs, blk) { }
    LU_piv_(Int_ type, 
            Int_ mt, Int_ mb, 
            Int_ rhs, 
            Int_ blk)  
      : LU_nopiv_<A_,B_>(type, mt, mb, rhs, blk) { }
    
    LU_piv_(const LU_piv_ &b) : LU_nopiv_<A_,B_>(b) { }

    using LU_nopiv_<A_,B_>::get_leading_flop;
    using LU_nopiv_<A_,B_>::decompose;
    using LU_nopiv_<A_,B_>::solve_1_x;
    using LU_nopiv_<A_,B_>::solve_2_x;
    using LU_nopiv_<A_,B_>::solve_2_r;
    
    virtual Int_ update() {

      if (this->get_mt() && this->get_mb()) {

#pragma omp task
        {
          apply_pivots( Mat_Base_::LEFT, Mat_Base_::NO_TRANSPOSE, 
                        this->_P, this->_ATR);
          trsm( Mat_Base_::LEFT, Mat_Base_::LOWER_TRIANGULAR,
                Mat_Base_::NO_TRANSPOSE, Mat_Base_::UNIT_DIAG,
                1, this->_ATL, this->_ATR );
        }
        
#pragma omp task
        {
          trsm( Mat_Base_::RIGHT, Mat_Base_::UPPER_TRIANGULAR,
                Mat_Base_::NO_TRANSPOSE, Mat_Base_::NONUNIT_DIAG, 
                1, this->_ATL, this->_ABL );
        }

#pragma omp taskwait
        
        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ABL, this->_ATR, 1, this->_ABR );

      }
      return true;
    }
    
    virtual Int_ decompose() {
      if (this->get_mt()) 
        lu_piv( this->_ATL, this->_P );

      this->update();

      return true;
    }

    virtual Int_ solve_1_x() {
      if (this->get_mt()) {
        apply_pivots( Mat_Base_::LEFT, Mat_Base_::NO_TRANSPOSE,
                      this->_P, this->_XT);
        trsm( Mat_Base_::LEFT, Mat_Base_::LOWER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::UNIT_DIAG,
              1, this->_ATL, this->_XT );
      }
      if (this->get_mt() && this->get_mb()) {
        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ABL, this->_XT, 1, this->_XB );
      }
      return true;
    }

    virtual Int_ solve_1_r() {
      if (this->get_mt()) {
        apply_pivots( Mat_Base_::LEFT, Mat_Base_::NO_TRANSPOSE,
                      this->_P, this->_RT);
        trsm( Mat_Base_::LEFT, Mat_Base_::LOWER_TRIANGULAR,
              Mat_Base_::NO_TRANSPOSE, Mat_Base_::UNIT_DIAG,
              1, this->_ATL, this->_RT );
      }
      if (this->get_mt() && this->get_mb()) {
        gemm( Mat_Base_::NO_TRANSPOSE, Mat_Base_::NO_TRANSPOSE,
              -1, this->_ABL, this->_RT, 1, this->_RB );
      }
      return true;
    }
  };

}

#endif
